

library(ggplot2)
library(tidyr)
setwd("C:/Users/at449/Desktop/replication package/matlab code")

#load the data
dt<-read.csv("data.csv", head = FALSE)
hist(dt$V1)


#Calculate the distribution of interest rate and LTV for simualtion
sorted_data <- dt[order(dt[, 1]), ]
totalBR<-sum(sorted_data$V5)
PrLTV75<-sum(sorted_data$V5[sorted_data$V1<=77.5])/totalBR
PrLTV80<-sum(sorted_data$V5[sorted_data$V1>77.5 & sorted_data$V1<=82.5 ])/totalBR
PrLTV85<-sum(sorted_data$V5[sorted_data$V1>82.5 & sorted_data$V1<=87.5 ])/totalBR
PrLTV90<-sum(sorted_data$V5[sorted_data$V1>87.5 & sorted_data$V1<=92.5 ])/totalBR
PrLTV95<-sum(sorted_data$V5[sorted_data$V1>92.5 & sorted_data$V1<=97.5 ])/totalBR

RLTV75<-sum(sorted_data$V4[sorted_data$V1<=77.5]*sorted_data$V5[sorted_data$V1<=77.5])/sum(sorted_data$V5[sorted_data$V1<=77.5])
RLTV80<-sum(sorted_data$V4[sorted_data$V1>77.5 & sorted_data$V1<=82.5 ]*sorted_data$V5[sorted_data$V1>77.5 & sorted_data$V1<=82.5 ])/sum(sorted_data$V5[sorted_data$V1>77.5 & sorted_data$V1<=82.5 ])
RLTV85<-sum(sorted_data$V4[sorted_data$V1>82.5 & sorted_data$V1<=87.5 ]*sorted_data$V5[sorted_data$V1>82.5 & sorted_data$V1<=87.5 ])/sum(sorted_data$V5[sorted_data$V1>82.5 & sorted_data$V1<=87.5 ])
RLTV90<-sum(sorted_data$V4[sorted_data$V1>87.5 & sorted_data$V1<=92.5 ]*sorted_data$V5[sorted_data$V1>87.5 & sorted_data$V1<=92.5 ])/sum(sorted_data$V5[sorted_data$V1>87.5 & sorted_data$V1<=92.5 ])
RLTV95<-sum(sorted_data$V4[sorted_data$V1>92.5 & sorted_data$V1<=97.5 ]*sorted_data$V5[sorted_data$V1>92.5 & sorted_data$V1<=97.5 ])/sum(sorted_data$V5[sorted_data$V1>92.5 & sorted_data$V1<=97.5 ])

distributionLTV<-c(PrLTV75,PrLTV80,PrLTV85,PrLTV90,PrLTV95)
distributionR<-c(RLTV75,RLTV80,RLTV85,RLTV90,RLTV95)

#Calculate the distribution of interest rate and LTV for other simualtion
#load the data
dt2<-read.csv("dataFC.csv", head = FALSE)
hist(dt2$V1)

FtotalBR<-sum(dt2$V5)
FPrLTV75<-sum(dt2$V5[dt2$V1<=77.5])/FtotalBR
FPrLTV80<-sum(dt2$V5[dt2$V1>77.5 & dt2$V1<=82.5 ])/FtotalBR
FPrLTV85<-sum(dt2$V5[dt2$V1>82.5 & dt2$V1<=87.5 ])/FtotalBR
FPrLTV90<-sum(dt2$V5[dt2$V1>87.5 & dt2$V1<=92.5 ])/FtotalBR
FPrLTV95<-sum(dt2$V5[dt2$V1>92.5 & dt2$V1<=97.5 ])/FtotalBR

FRLTV75<-sum(dt2$V4[dt2$V1<=77.5]*dt2$V5[dt2$V1<=77.5])/sum(dt2$V5[dt2$V1<=77.5])
FRLTV80<-sum(dt2$V4[dt2$V1>77.5 & dt2$V1<=82.5 ]*dt2$V5[dt2$V1>77.5 & dt2$V1<=82.5 ])/sum(dt2$V5[dt2$V1>77.5 & dt2$V1<=82.5 ])
FRLTV85<-sum(dt2$V4[dt2$V1>82.5 & dt2$V1<=87.5 ]*dt2$V5[dt2$V1>82.5 & dt2$V1<=87.5 ])/sum(dt2$V5[dt2$V1>82.5 & dt2$V1<=87.5 ])
FRLTV90<-sum(dt2$V4[dt2$V1>87.5 & dt2$V1<=92.5 ]*dt2$V5[dt2$V1>87.5 & dt2$V1<=92.5 ])/sum(dt2$V5[dt2$V1>87.5 & dt2$V1<=92.5 ])
FRLTV95<-sum(dt2$V4[dt2$V1>92.5 & dt2$V1<=97.5 ]*dt2$V5[dt2$V1>92.5 & dt2$V1<=97.5 ])/sum(dt2$V5[dt2$V1>92.5 & dt2$V1<=97.5 ])

distributionLTVF<-c(FPrLTV75,FPrLTV80,FPrLTV85,FPrLTV90,FPrLTV95)
distributionRF<-c(FRLTV75,FRLTV80,FRLTV85,FRLTV90,FRLTV95)



#Calculate the distribution of interest rate and LTV for data
dtI<-read.csv("dataI.csv", head = FALSE)
dtI$V5<-1
ItotalBR<-sum(dtI$V5)
IPrLTV75<-sum(dtI$V5[dtI$V1<=77.5])/ItotalBR
IPrLTV80<-sum(dtI$V5[dtI$V1>77.5 & dtI$V1<=82.5 ])/ItotalBR
IPrLTV85<-sum(dtI$V5[dtI$V1>82.5 & dtI$V1<=87.5 ])/ItotalBR
IPrLTV90<-sum(dtI$V5[dtI$V1>87.5 & dtI$V1<=92.5 ])/ItotalBR
IPrLTV95<-sum(dtI$V5[dtI$V1>92.5 & dtI$V1<=97.5 ])/ItotalBR


IRLTV75<-mean(dtI$V4[dtI$V1<=77.5])
IRLTV80<-mean(dtI$V4[dtI$V1>77.5 & dtI$V1<=82.5 ])
IRLTV85<-mean(dtI$V4[dtI$V1>82.5 & dtI$V1<=87.5 ])
IRLTV90<-mean(dtI$V4[dtI$V1>87.5 & dtI$V1<=92.5 ])
IRLTV95<-mean(dtI$V4[dtI$V1>92.5 & dtI$V1<=97.5 ])

distributionLTVI<-c(IPrLTV75,IPrLTV80,IPrLTV85,IPrLTV90,IPrLTV95)
distributionRI<-c(IRLTV75,IRLTV80,IRLTV85,IRLTV90,IRLTV95)

hist(dtI$V1)

plot(distributionLTV)



####### Plot distribution of the LTV

test_data <-
  data.frame(
    SocialPlanner = distributionLTV,
    InitialData = distributionLTVI,
    LTV = c(75,80,85,90,95)
  )
test_data %>%
  gather(key,value, SocialPlanner, InitialData) %>%
  ggplot(aes(x=LTV, y=value, colour=key)) +
  geom_point()+ ggtitle("Market share by LTV")

ggsave("MarketShareLTV.png")


####### Plot distribution of the Interst rate

test_data <-
  data.frame(
    SocialPlanner = distributionRI,
    InitialData = distributionR,
    LTV = c(75,80,85,90,95)
  )
test_data %>%
  gather(key,value, SocialPlanner, InitialData) %>%
  ggplot(aes(x=LTV, y=value, colour=key)) +
  geom_point()+ ggtitle("Average interest rate by LTV")

ggsave("RateLTV.png")

################################

####### Plot distribution of the LTV

test_data <-
  data.frame(
    SocialPlanner = distributionLTVF,
    InitialData = distributionLTVI,
    LTV = c(75,80,85,90,95)
  )
test_data %>%
  gather(key,value, SocialPlanner, InitialData) %>%
  ggplot(aes(x=LTV, y=value, colour=key)) +
  geom_point()+ ggtitle("Market share by LTV")

ggsave("FMarketShareLTV.png")


####### Plot distribution of the Interst rate

test_data <-
  data.frame(
    SocialPlanner = distributionRI,
    InitialData = distributionRF,
    LTV = c(75,80,85,90,95)
  )
test_data %>%
  gather(key,value, SocialPlanner, InitialData) %>%
  ggplot(aes(x=LTV, y=value, colour=key)) +
  geom_point()+ ggtitle("Average interest rate by LTV")

ggsave("FRateLTV.png")





